*Figure 1 - factor structure of attitudes toward divine

factor Absolute-Never_fails, pf
rotate, oblimin blanks(.5)
loadingplot, ylabel(0 "0.0" 0.2 "0.2" 0.4 "0.4" 0.6 "0.6" 0.8 "0.8" 1.0 "1.0", labcolor(black) grid glcolor(black%20) glpattern(solid) nogextend) ///
 xlabel(0 "0.0" 0.2 "0.2" 0.4 "0.4" 0.6 "0.6" 0.8 "0.8" 1.0 "1.0", labcolor(black) grid glcolor(black%20) glpattern(solid) nogextend) ///
 plotregion(fcolor(white)) ///
 xscale(r(-.05 1.05) lc(black)) yscale(r(-.05 1.05) lc(black)) ///
 ytitle("Negative qualities") ///
 xtitle("Positive qualities")
 
*Figure 3 - comparing LCA solution with factor analysis

twoway(scatter LCA pos01, mcolor(white) mlcolor(black) ylabel(1 `""Class 1: Indefinite" "and unengaged""' 2 `""Class 2: Loving," "not judgmental," "and engaged""' 3 `""Class 3: God" "as Nullity""' 4 `""Class 4: Loving," "not judgmental," "and unengaged""' 5 `""Class 5: Loving," "judgmental," "and engaged""')) || scatter latent_divine neg01, ///
 mcolor(black) ///
 ylabel(,labcolor(black) grid glcolor(black%20) glpattern(solid)) ///
 xlabel(0 "0.0" 0.2 "0.2" 0.4 "0.4" 0.6 "0.6" 0.8 "0.8" 1.0 "1.0",labcolor(black) grid glcolor(black%20) glpattern(solid)) ///
 xtitle("Factor scores where higher values = stronger association with those qualities") ///
 plotregion(fcolor(white)) ///
 xscale(r() lc(black)) yscale(r() lc(black)) ///
 
*Figure 4 - composition of classes by trad denom

histogram denom, by(latent_divine) percent discrete ///
 ylabel(0 20 40 60,labcolor(black) grid glcolor(black%20) glpattern(solid) nogextend) ///
 xlabel(1 "Mainline Protestant" 2 "Catholic" 3 "Jewish" 4 "Other" 5 "None" 6 "Evangelical" 7 "Black Protestant",labcolor(black) grid glcolor(black%20) glpattern(solid)) ///
 plotregion(fcolor(white)) ///
 yscale(r() lc(black))
 
*Figure 5 - expressive religious behaviors / literalism

mean performance01, over(latent_divine)
est store perf
mean literalism01 , over(latent_divine)
est store literal

coefplot perf literal, ylabel(1 `""Class 1: Indefinite" "and unengaged""' 2 `""Class 2: Loving," "not judgmental," "and engaged""' 3 `""Class 3: God" "as Nullity""' 4 `""Class 4: Loving," "not judgmental," "and unengaged""' 5 `""Class 5: Loving," "judgmental," "and engaged""', labcolor(black) grid glcolor(black%20) glpattern(solid) nogextend) xlabel(,labcolor(black) grid glcolor(black%20) glpattern(solid) nogextend)  plotregion(color(white)) ///
 levels(95) ///
 xscale(r() lc(black)) yscale( lc(black))

 *Figure 6 - predicted marginal effects of covariates on class inclusion
 
mlogit latent_divine evangel blackprot none jewish catholic other charismatic performance01 literalism01 universalism01 trustscale01 control01  age white black educ male income[pweight=weight], base(4)
margins, dydx(*) predict(outcome(1)) post
est store all1
mlogit latent_divine evangel blackprot none jewish catholic other charismatic performance01 literalism01 universalism01 trustscale01 control01  age white black educ male income[pweight=weight], base(4)
margins, dydx(*) predict(outcome(2)) post
est store all2
mlogit latent_divine evangel blackprot none jewish catholic other charismatic performance01 literalism01 universalism01 trustscale01 control01  age white black educ male income[pweight=weight], base(4)
margins, dydx(*) predict(outcome(3)) post
est store all3
mlogit latent_divine evangel blackprot none jewish catholic other charismatic performance01 literalism01 universalism01 trustscale01 control01  age white black educ male income[pweight=weight], base(4)
margins, dydx(*) predict(outcome(4)) post
est store all4
mlogit latent_divine evangel blackprot none jewish catholic other charismatic performance01 literalism01 universalism01 trustscale01 control01  age white black educ male income[pweight=weight], base(4)
margins, dydx(*) predict(outcome(5)) post
est store all5

coefplot all1, bylabel("Class 1: Indefinite," "and unengaged") || all2, bylabel("Class 2: Loving," "not judgmental," "and engaged") || all3, bylabel("Class 3: God" "as Nullity" " ") || all4, bylabel("Class 4: Loving," "not judgmental," "and unengaged") || all5, bylabel("Class 5: Loving," "judgmental," "and engaged") ||, xline(0) groups(evangel other = "{bf: Denomination}" charismatic universalism01 = "{bf:Relig belief, practice}" trustscale01 income = "{bf: Demographics}") levels(95 90) ///
plotregion(fcolor(white) lcolor(black) lp(solid))  byopts(row(1)) coeflabel(evangel="Evangelical" blackprot="Black protestant" none="None" jewish="Jewish" catholic="Catholic" other="Other relig" charismatic="Charismatic" performance01="Performance" literalism01="Bib literalism" universalism01="Universalism" trustscale01="Social trust" control01="Perceived control" educ="Education" age="Age" white="White" black="Black" male="Male" income="Income") ///
ylabel(,labcolor(black) grid glcolor(black%20) glpattern(solid)) norecycle legend(off) nokey ///
xlabel(-1 -.5 0 .5 ,labcolor(black) grid glcolor(black%20) glpattern(solid)) 

